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FIGURE  CAPTIONS 


Figure  1.  Semtdiuraal  density  amplitude  vs.  height  and  ladtiule  at  42*N. 

Figure  2.  Percent  density  variation  vs.  height  and  local  time  over  the  equator 
during  March. 

Figure  3.  Semidiurnal  (top)  and  diurnal  (bottom)  density  variations  as  a  function 
of  height  near  the  equator,  for  the  1979  Kwajalein  reference  atmosphere  (solid  lines) 
and  densities  from  the  DCORR  subroutine. 

Figure  4.  Daytime  densities  at  200  km  as  a  result  of  one-unit  Kp  binning  and  10* 
latitude  binning  over  the  20-day  observational  period.  Top:  SETA  densities 
(normalized  to  200  km)  versus  Kp  at  various  geomagnetic  latitudes  (left) 
and  versus  geomagnetic  latitude  for  Kp  =  4.0  and  6.0  (right). 

Middle:  Same  as  top  figures,  except  for  the  TIGCM  at  200  km. 

Bottom:  Same  as  above,  except  for  the  MSISE90  model  at  200  km. 

Figure  5.  Predictions  for  the  storm  response  in  TEC  at  Hamilton,  MA  for  the 
23  November,  1987  storm,  compared  with  the  observations  (solid  line). 

Figure  6.  (a)  Simulation  results  for  the  event  on  4  February  1984  at  Ramey, 
together  with  observations,  (b)  The  required  vertical  drifts  (V^,  positive  for  upward) 
at  350  km,  induced  by  a  disturbance  electric  field  (E^  positive  for  eastward). 

The  dotted  vertical  lines  are  used  to  indicate  the  temporal  correlation  between 
the  ionospheric  patterns  and  the  enhancement  and  recovery 
phases  of  substorms. 

Figure  7.  The  mean  diurnal  variations  of  TEC  for  all  12  months  of  1983 
at  Ramey,  Puerto  Rico  (solid  line),  compared  to  the  predictions  from  models: 
LMTEC  (dashed),  Damon-Hartranft  (dot-daslred).  Bent  (dot-dot-dot-dashed) 
and  IRI-90  (long-dashed). 

Figure  8.  Examples  of  fits  to  Millstone  Hill  IS  radar  profiles  using  a  double 
variable  Chapman  layer  formalism. 


PREFACE 


Contract  F19628-90-K-0014  consisted  of  two  main  task  areas,  NEUTRAL  DENSITY  MODEL¬ 
ING  and  PLASMA  DENSITY  MODELING.  Work  completed  under  these  categories  are  reported 
below  in  Sections  1  and  2,  respectively. 


1.  NEUTRAL  DENSITY  MODELING 


1.1  Technique  for  Assimilative  Data  Anidysis 

In  this  effort  a  set  of  tabulated  functions  called  "Hough  Mode  Extensions"  (HMEs),  which  repre¬ 
sent  numerical  extensions  of  classical  Hough  modes  into  the  viscous  regime  of  the  thermosphere, 
are  used  to  least-squares  fit  a  climatological  data  base  of  tidal  measurements.  The  data  base  con¬ 
sists  of  monthly  average  vertical  profiles  of  semidiurnal  amplitudes  and  phases  at  seventeen  radar 
sites  accessing  some  part  of  the  80  - 1  SO  km  height  region.  The  radars  are  distributed  between  78*S 
and  70**N  latitude,  and  each  one  provides  measurements  of  one  or  more  of  the  following:  eastward 
wind,  southward  wind,  perturbation  temperature.  As  a  result  of  the  fitting  process,  a  single  com¬ 
plex  ntHmalizing  coefficient  is  derived  for  each  month  and  for  each  of  the  four  HMEs,  (tesignated 
(2,2),  (2,3),  (2,4)  and  (2,5)  after  their  classical  Hough  function  designations.  Once  the  complex 
coefficients  are  derived,  reconstruction  by  weighted  superposition  of  the  HMEs  results  in  globally 
continuous  specifications  of  semidiurnal  horizontal  and  vertical  wind,  temperature,  pressure,  and 
density  throughout  the  80  -  150  km  height  region.  The  tidal  variations  in  density,  in  particular, 
provide  greater  accuracy  for  several  aerospace  applications.  An  example  of  the  semidiurnal  den¬ 
sity  amplitudes  obtained  by  this  method  are  depicted  in  Figure  1.  A  similar  data  calibration  to 
numerically-generated  tidal  functions  was  performed  for  the  diurnal  tide.  When  combined  with  the 
semidiurnal  tide,  one  is  able  to  obtain,  for  any  month  and  any  latitude,  the  type  of  percent  density 
variation  plots  illustrated  in  Figure  2  for  0*  latitude  during  the  month  of  March. 

The  methodology  described  here  can  also  be  utilized  to  derive  tidal  lower  boundary  conditions 
for  Thermospheric  General  Circulation  Models  (TGCMs),  or  as  a  basis  for  future  empirical  moctel 
development 

A  complete  description  of  the  above  work  is  contained  in  the  manuscript  entitled 

SEMIDIURNAL  TIDE  IN  THE  80  - 130  KM  REGION:  AN  ASSIMILATIVE  DATA  ANAL¬ 
YSIS,  by  J.M.  Forbes,  A.H.  Manson,  R.A.  \%cent,  G.J.  Fraser,  S.K.  Avery,  R.R.  Clark,  R.  Sch- 
minder,  and  D.  Kurschner. 

which  has  been  accepted  for  publication  in  the  Journal  of  Atmospheric  and  Terrestrial  Physics,  a 
preprint  copy  of  which  has  been  delivered  to  PL/GPIM.  In  this  publication,  comparisons  are  also 
made  with  HME  coefficients  and  global  tidal  fields  from  the  Forbes  and  Vial  [  1989]  numerical  tidal 
model,  and  with  independent  data  sets  not  included  in  the  least-squares  fitting  process. 

12  Density  Correction  Algorithm  for  Groves/MSIS  Model 
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As  a  result  of  the  fitting  procedure  described  jd)ove,  the  diurnal  and  semidiurnal  perturbation  densi¬ 
ties  have  been  fit  analytically  with  Legendre  polynomials,  and  an  algorithm  (consisting  of  basically 
three  FORTRAN  subroutines)  developed  which  uses  these  expansions  to  provide  a  density  correc¬ 
tion  to  the  so-called  ”Groves/MSIS83’'  density  model  in  operational  use  by  the  Air  Weather  Service. 
A  brief  description  of  the  GMSIS83  test  program,  aiul  the  TMOD.  DCORR.  and  DHME  subrou¬ 
tines  developed  under  this  effort  are  described  below,  and  are  provided  in  full  in  the  APPENDIX. 

GMSIS83 

This  is  a  demonstration  program  which  calculates  the  percent  mass  density  variation  about  diur¬ 
nal  mean  values  at  various  heights  and  latitudes  using  (a)  the  Groves/MSIS83  Model  and  (b)  the 
DCORR  subroutine  developed  under  the  present  contract  Tables  of  percent  deviation  vs.  height 
(80  -  200  km)  and  local  time  (0  -  24  hours)  are  generated  at  latitudes  of  0. 30.  and  60  degrees.  The 
intent  is  to  demonstrate  use  of  the  DCORR  subroutine  to  those  familiar  with  use  of  the  GROVES 
subroutine  (i.e..  the  one  currently  in  use  by  AWS). 

In  the  GMSIS83  code,  density  variations  are  first  calculated  with  GROVES.  Further  below,  a 
similar  code  appears,  except  that  the  percent  total  mass  density  variation  is  provided  through 
a  call  to  TMOD  (see  below)  which  subsequently  calls  DCORR.  For  altitudes  between  145  and 
165  km.  TMOD  provides  a  linear  interpolation  between  the  DCORR  results  at  145  km.  and  the 
Groves/MSIS83  results  at  165  km.  so  that  no  density  discontinuity  exists  as  a  result  of  using 
DCORR. 

It  is  important  to  note  that  DCORR  only  corrects  the  tidal  density  variation;  that  is.  when  it  is  imple¬ 
mented  in  the  illustrated  fashion,  as  a  correction  to  the  diurnal  mean  density  from  Groves/MSIS83. 
the  other  parameters  (temperature,  pressure,  etc.)  given  by  Groves/MSlS83  are  the  diurnal  mean 
values  with  no  tidal  variation  whatsoever.  If  such  variations  are  required,  the  Groves/MSIS83 
model  would  have  to  be  called  again  with  the  appropriate  input  parameters. 

TMOD 


This  subroutine  provides  Groves/MSIS83  densities  with  "corrected”  tidal  perturbations  between 
80  and  145  km.  including  a  smooth  merger  with  Groves/MSIS83  between  145  and  165  km. 

DCORR 

Subroutine  DCORR  computes  the  density  correction  factor  dens  =  dens*rc  to  account  for  density 
perturbations  due  to  tides  in  the  80  to  150  km  height  regime.  The  correction  factor  is  determined 
from  Legendre  polynomial  reconstruction  of  diurnal  (n=l)  and  semidiurnal  (n=2)  harmonics  of 
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perturbation  density.  This  subroutine  also  provides  relative  density  amplitude  {Ap/p)  and  phase 
(local  time  of  maximum)  for  the  n=l  and  n=2  harmonics. 

This  subroutine  contains  the  normalizing  coefficients  for  the  Hough  Mode  Extension  (HME)  func¬ 
tions  (see  Section  1.1),  obtained  by  calls  to  the  DHMEnm  (where  n=l,  m=l  and  n=2,  m=2,3,4,5) 
subroutines.  The  coefficients  (one  set  for  each  month  of  the  year)  are  |HX)vided  in  data  statements. 
Two  options  are  provided  for  the  semidiurnal  tide:  (a)  from  the  Fm'bes  and  [1989]  numeri¬ 
cal  model;  and  (b)  from  least-squares  fitting  of  a  global  array  of  wind  and  temperature  data  from 
ground-based  radars  (see  Section  1.1).  The  latter  set,  being  derived  from  climatlogical  averages 
of  actual  measurements,  is  the  default  in  the  present  version  of  the  code.  The  table  of  coefficients 
derived  from  data  should  be  updated  as  more  measurements  and  validation  results  become  avail¬ 
able.  Operationally,  the  tidal  density  correction  algorithm  should  be  implemented  so  that  these 
coefficients  can  easily  be  updated. 

DHMEnm  where  n=l,  m=l  or  ns2,  ms23»4,S 

These  subroutines  interpolate  tables  of  6-th  order  Legendre  polynomial  coefficients  to  give  the 
complex  HME  for  perturbation  density  (percent  of  mean)  of  the  (n.m)  mode. 

A  3.5-inch  DOS  compatible  disk  containing  the  files  GMSIS83,  GROVES,  DCORR,  and  TMOD 
were  sent  to  Captain  Mark  Raffensberger  at  HQ  Air  Weather  Service/XTX,  and  to  Mr.  EA.  Marcos, 
technical  monitor  of  this  contract  Documentation  and  other  descriptive  information  regarding  the 
methodology  of  constructing  the  algorithm  were  also  sent. 

The  availability  of  density  data  distributed  over  sufficient  local  times  to  allow  tidal  determina¬ 
tions  of  density  within  the  80  -  150  km  altitude  regime,  and  true  validation  of  the  DCORR  algo¬ 
rithm,  are  extremely  scarce.  One  data  source  can  be  found  in  "Kwajalein  Reference  Atmosf^res, 
1979”,  AFGL-TR-79-0241, 24  September,  1979.  Here  analyses  of  seventeen  high-altitude  ROBIN 
spheres,  launched  within  48  hours  on  19  to  21  July  1978  at  Kwajalein,  provide  data  for  determina¬ 
tion  of  diurnal  and  semidiurnal  variations  in  density.  The  corresponding  diurnal  and  semidiurnal 
amplitudes  are  plotted  vs  height  as  solid  lines  in  Figure  3,  with  the  DCORR  output  indicated  by 
dots.  The  diurnal  amplitudes  agree  reasonably  with  the  Kwajalein  values,  but  the  semidiurnal  am¬ 
plitudes  are  about  half  of  the  Kwajalein  values.  However,  one  should  not  put  too  much  weight  into 
the  above  comparison  as  it  suffers  from  a  number  of  shortcomings.  The  Kwajalein  data  set  may 
not  be  very  representative,  and  uncertainties  exist  due  to  measurement  technique  and  to  the  fourier 
data  analysis;  however,  no  imformation  on  these  possible  sources  of  error  are  given  in  the  report 
On  the  other  hand,  real  differences  may  exist  at  Kwajalein  due  to  longitude  effects  which  are  not 
taken  into  account  in  the  DCORR  formulation.  A  meaningful  evaluation  of  DCORR  requires  data 
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that  are  not  yet  in  existence.  However,  the  thetMetical  basis  of  the  present  method  is  reasonably 
sound,  and  offers  the  best  possible  means  of  translating  a  large  data  base  of  wind  measurements 
into  density  variations  which  are  not  piesendy  possible  to  measure  on  a  global  scale. 

13  SETA  DensityAVind  Analysis  and  Model  Validations 

Under  this  effort,  satellite-based  measurements  were  utilized  to  elucidate  the  latitudinal,  local¬ 
time,  and  magnetic  activity  dependence  of  winds  and  densities  in  the  scantily-observed  atmospheric 
region  between  170  and  220  km  above  45*  magnetic  latitude.  One  data  set  consists  of  atmosi^ric 
densities  from  high-accuracy  (time  resolution  »  5-6  hours)  orbital  analyses  of  thr^  Doppler  Beacon 
satellites  in  orbit  during  1973.  The  perigees  of  these  satellites  are  generally  restricted  to  160-180 
km,  1200-1400  LST,  and  geographic  latitudes  greater  than  about  45*.  Statistical  relationships  are 
derived  between  the  density  changes  and  the  planetary  magnetic  index.  Kp,  and  the  five-hour  mean 
of  the  auroral  electrojet  index.  XE.  The  former  relationships  are  compared  with  those  derived  from 
the  MSISE90  empirical  model  [Hedin,  1991],  which  is  found  to  overestimate  the  rate  of  increase 
of  density  with  respect  to  Kp. 

Another  data  set  consists  of  densities  and  cross-track  winds  from  the  Satellite  Elecu-ostatic  Triax- 
ial  Accelerometer  (SETA)  experiment  for  the  21  March  -  9  April,  1979,  period,  which  includes 
several  intervals  of  elevated  magnetic  activity.  Besides  comparing  various  time  series,  the  data 
are  also  binned  according  to  10*  latitude  increments  and  unit  increments  of  Kp  to  derive  trends. 
Some  typical  results  include  the  following,  corresponding  to  average  changes  in  the  45  to  65*  mag¬ 
netic  latitude  band  as  Kp  is  increased  from  1  to  6:  (1)  For  the  nightside  (»  2230  LT),  a  change  in 
cross-track  (nearly  zonal)  wind  from  25  ±  25  ms“*  (eastward)  to  -125  ±  25  ms"*  (westward),  and 
an  increase  of  about  20%  in  density;  (2)  For  dayside  (»  1030  LT),  a  change  in  cross-track  wind 
from  25  ±  25  ms”*  (eastward)  to  125  ±  25  ms"*  (eastward),  and  a  density  increase  of  25%.  For 
some  individual  sudden  enhancements  in  magnetic  activity,  changes  in  winds  and  densities  can  be 
more  than  double  the  above  average  values.  Comparisons  are  also  made  with  the  NCAR  TIGCM 
(National  Center  for  Atmospheric  Research  Thermosphere-Ionosphere  General  Circulation  Model) 
simulation  for  the  complete  20-day  interval,  and  with  recent  empirical  models  of  densities  [Hedin, 
1991]  and  winds  [Hedin  et  al.,  1991],  with  data  points  in  each  case  derived  for  the  satellite  paths, 
instrument  orientations,  and  sampling  rates  identical  to  that  of  the  SETA  experiment 

A  typical  example  of  the  type  of  comparisons  made  in  this  study  are  provided  in  Figure  4.  Figure 
4  illustrates  daytime  («  1030  LT)  density  vs.  Kp  fcnr  several  latitudes  (left  panels)  for  SETA  (top), 
TIGCM  (middle),  and  MSISE90  (bottom),  and  also  the  latitudinal  dependences  of  density  (right 
panels)  for  the  Kp  =  4  and  Kp  =  6  bins.  No  time  delay  is  assumed  in  the  SETA  analysis  and 
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comparisons.  Note  that  the  absolute  scales  for  the  TIGCM  densities  are  slightly  lower,  due  to  the 
net  DC  level  deficiency  of  the  TIGCM.  However,  as  indicated  by  the  slopes  of  the  curves  in  the 
left-hand  panels,  the  relative  increase  of  density  with  Kp  is  very  nearly  the  same  ftn*  the  all  three 
models.  The  slopes  for  SETA  and  TIGCM  are  more  or  less  independent  of  latitude  up  to  Kp  = 
4.  However,  for  higher  Kp  the  density  increase  (slope)  is  greatest  at  70”  and  least  at  80*.  for  both 
the  TIGCM  and  SETA  analyses.  This  shows  up  as  maxima  in  the  density  vs.  latitude  plots  (right 
panels)  near  70*  magnetic  latitude,  reflecting  the  maximum  in  joule  heating  occurring  in  the  daytime 
auroral  oval.  Adiabatic  heating  and  cooling  effects  connected  with  the  dynamical  response  of  the 
thermosphere  represent  alternative  processes  which  can  give  rise  to  such  maxima  (or  minima)  in 
the  density  field  [Crowley  et  al.,  1989;  Schoendorf  et  al..  1991].  The  MSISE90  model  suggests 
only  a  slight  tendency  for  maximum  density  response  at  70*  latitude  for  Kp  ^  6.  as  indicated  in  the 
lower  left  panel  of  Figure  4. 

A  complete  description  of  the  work  completed  under  this  effort  is  contained  in  the  mansucript; 

MAGNETIC  ACTIVITY  DEPENDENCE  OF  HIGH-LATITUDE  THERMOSPHERIC 
WINDS  AND  DENSITIES  BELOW  200  KM.  by  J.  M.  Forbes,  R.  G.  Roble,  and  F.  A  Marcos. 

which  has  been  accepted  for  publication  in  the  Journal  of  Geophysical  research,  a  preprint  copy  of 
which  has  been  delivered  to  PL/GPIM. 

2.  PLASMA  DENSITY  MODELING 

Ongoing  developments  and  advances  in  the  held  of  ionsopheric  modeling  come  from  a  variety  of 
sources  and  motivations.  Some  of  these  are  in  response  to  the  concerns  of  the  user-community, 
whereby  convenient  descriptions  of  ionospheric  phenomena  are  required  so  that  systems  can  be 
adjusted  to  account  for  both  long-term  average  behavior  and  short-term  perturbations.  The  main 
issue  addressed  in  the  research  presented  here  is 

-  can  better  descriptions  be  found  ? 

A  number  of  different  analyses  are  presented  herein.  Broadly  speaking  these  fit  into  two  categories: 

Data-based  This  type  of  study  involves  analysis  of  new  or  revised  data  sets,  to  determine  new 
trends  in  the  data,  or  to  test  and  validate  a  new  numerical  description  of  data,  providing  an  improved 
description  of  variations  in  the  data  to  the  user. 
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Modd'based  This  type  of  study  involves  the  use  of  a  theoretical  ionospheric  model  used  as  a 
diagnostic  tool  to  study  particular  ionospheric  phenomena.  This  serves  both  as  a  validation  step  on 
the  theoretical  model,  and  at  the  same  times  imj^oves  our  physical  understanding  of  the  processes 
involved,  providing  in  atm  a  greater  level  of  “predictability"  to  the  user. 

2.1  Day«t(Mlay  Variability  and  Stonn«time  Response 

In  this  effort,  the  emphasis  was  to  quantify,  and  where  possible,  describe,  the  day  to  day  variability 
observed  in  the  ionosphere. 

Two  solar  cycles  worth  of  hourly  data  from  both  Hamilton.  MA  (total  content,  or  TEC)  and  Wallops 
Is.,  VA  (Nmax)  were  used  as  the  basis.  Slab  thickness  values  were  evaluated  making  only  a  small 
correction  in  local  time  between  the  two  locations.  Firstly,  the  degree  of  day-to-day  variability 
was  quantified  on  a  monthly  basis,  and  was  seen  to  be  as  significant  as  the  degree  of  agreement 
currently  achievable  with  existing  climatological  models.  Secondly,  the  level  of  predictability  of 
these  variations  were  studied,  employing  the  correlation  between  pairs  of  these  parameters  and 
testing  the  validity  of  “persistence”  forecasting  over  various  periods.  Overall,  the  bulk  of  the  vari¬ 
ations  were  not  accessible  by  these  simple  means.  Next,  it  was  established  that  the  periods  of 
greatest  variability  were  shown  to  be  geomagnetically  active  periods,  and  a  separate  study  of  the 
400  geomagnetic  storms  in  the  database  was  undertaken.  The  average  storm-time  variations  were 
described  and  related  to  typical  physical  processes,  but  at  the  same  time  it  was  shown  that  a  typ¬ 
ical  storm  and  an  average  storm  are  quite  different.  A  second  model  of  ionospheric  reponse  in 
TEC  was  developed  based  on  probabilistic  considerations  of  individual  storm-time  features,  and 
this  achieved  quite  good  results,  qualitatively,  in  predicting  ionospheric  response.  An  example  of 
predicted  ionospheric  response  is  shown  in  Figure  5,  in  both  the  absolute  TEC  and  the  percentage 
deviation  from  expected  behavior. 

A  complete  description  of  this  work  is  contained  in  the  report 

THE  DAY-TO-DAY  VARIABILITY  OF  TOTAL  CONTENT,  PEAK  DENSITY  AND  SLAB 
THICKNESS,  AND  THE  IONOSPHERIC  RESPONSE  TO  GEOMAGNETIC  STORMS,  by 
M.W.  Fox  (1990) 

that  has  been  presented  to  PL/GPIM. 

2.2  A  New  Global  Ionospheric  Specification  Model 

It  has  long  been  a  goal  of  ionospheric  researchers  to  provide  a  realistic  description  of  the  global 
ionosphere  that  is  accessible  to  the  user  community  and  responsive  to  real-time  data.  Empirical 
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models  of  various  types  (e.g.  Llewellyn  and  Bent,  1973;  Rawer  and  Piggott,  1988)  provide  conve¬ 
nient  summaries  of  availiable  data  yet  rely  on  numerical  interpolations  in  locations  and  conditions 
not  covered  in  the  initial  database.  Theoretical  models  typically  involve  a  large  amount  of  time  to 
generate  an  ionosphere  for  a  given  set  of  conditions,  so  are  not  as  accessible. 

The  Parameterized  Real-time  Ionospheric  Specification  Model,  or  PRISM,  involving  a  significant 
collaboration  of  workers,  is  one  attempt  to  address  this  situation.  Theoretical  models  appropriate  to 
specific  geophysical  regions,  with  the  Phillips  Laboratory  model  (Anderson,  1973)  being  the  main 
component,  have  been  run  over  a  comprehensive  grid  of  conditions,  and  the  resultant  variations 
summarized  in  a  convenient  numerical  form.  Empirical  Orthonormal  Functions  (EOFs)  form  the 
basis  of  electron  density  profile  specification.  This  simple  parameterized  model  forms  the  basis 
for  ionospheric  description  in  PRISM.  The  second  stage  of  PRISM  involves  a  response  to  avail¬ 
able  data,  whereby  the  global  ionosphere  is  adjusted  according  to  the  departures  seen  from  the 
parameterized  background  model  and  the  observations. 

A  complete  description  of  this  work  and  of  the  model  validation  is  contained  in  the  report 

PRISM:  A  PARAMETERIZED  REAL-TIME  IONOSPHERIC  SPECIHCATION  MODEL 
by  R.E.  Daniell,  W.G.  Whaitenby,  L.D.  Brown,  D.N.  Anderson,  M.W.  Fox,  P.H.  Doherty,  D.T. 
Decker,  J.J.  Sojka  and  R.W.  Schunk  (1993) 

currently  under  preparation  at  PL/GPIM. 

23  A  Diagnostic  Modeling  Study  :  Diurnal  Twin  Peaks 

This  study  comprised  an  investigation  of  one  particular  type  of  diurnal  variation  in  total  electron 
content,  namely  the  diurnal  twin-peak.  Initially,  hourly  data  from  three  sites  in  the  Eastern  US 
sector,  provided  by  PL/GPIM,  was  analyzed  to  characterize  the  nature  of  this  phenomenon,  de¬ 
scribing  its  frequency,  magnitude  and  variations,  in  addition  to  the  correlation  between  the  feature 
at  the  three  locations.  It  was  then  established  that  the  feature  was  both  more  frequent  and  of  greater 
magnitude  during  periods  of  geomagnetic  distuitances,  with  the  twin  peaks  being  related  to  sub¬ 
storm  signatures.  Case  studies  suggested  that  the  likely  dynamical  source  for  this  phenomenon  are 
penetration  electric  fields,  as  demonstrated  by  the  time  scale  for  propagation  of  the  effect.  Theo¬ 
retical  modeling  studies  of  twin-peak  phenomena  in  turn  suggested  that  a  combination  of  latitude- 
(equivalently,  altitude- )dependent  electric  fields  and  chemical  loss  effects  are  responsible  for  the 
magnitude  of  the  effects  observed.  Figure  6  shows  an  example  of  the  modeling  results  and  the  na¬ 
ture  of  the  additional  vertical  drifts/electric  fields  required  to  match  the  observations  in  a  particular 
case. 
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A  complete  description  of  this  work  is  contained  in  the  report 

DIURNAL  DOUBLE  MAXIMA  IN  THE  F  REGION  IONOSPHERE:  SUBSTORM- 
RELATED  ASPECTS  by  X.Q.  Pi.  M.  MendUlo.  M.W.  Fox  and  D.N.  Anderson  (1993) 

that  has  been  delivered  to  PL/GPIM. 

2.4  Modeling  Ionospheric  Nightglow 

A  new  series  of  instruments  are  being  designed  and  built  for  future  usage  by  the  Defense  Meteo¬ 
rological  Satellite  Program,  DMSP.  These  will  include  ultraviolet  imaging  sensors  that  will  make 
observations  of  nightglow  both  at  the  nadir  and  at  the  earth’s  limb.  The  effwt  conducted  at  Boston 
University  consisted  of  applying  the  Phillips  Laboratory  ionospheric  model  (the  basis  of  the  PRISM 
model,  described  above)  at  low  and  middle  latitudes  to  generate  a  system  of  calibrations  between 
the  observed  nightglow  intensities  and  key  ionospheric  profile  paramaters  below  the  satellite,  thus 
enabling  ionospheric  monitoring  from  DMSP.  The  ionospheric  model  was  run  over  a  wide  grid  of 
conditions  (encompassing  the  grid  used  for  PRISM),  but  with  additional  variation  of  the  highly 
variable  vertical  E  x  B  drift  patterns  at  low  latitudes  in  order  to  build  into  the  system  some  ro¬ 
bustness  with  respect  with  day-to-day  variability.  This  model  ionosphere  was  then  combined  with 
the  neutral  atmosphere  from  MSIS-86  (Hedin,  1987),  and  the  “best  available”  set  of  reaction  rates 
that  describe  the  emissions  at  1356  X  (through  radiative  and  ion-ion  recombination)  and  6300  X 
(through  dissociative  recombination).  The  calibrations  were  then  derived  by  making  multiple  point- 
by-point  comparisons  between  the  derived  airglow  parameters  and  the  parameters  of  the  original 
ionosphere.  Sets  of  coefficients  have  been  established  that  relate  the  observed  nightglow  inten¬ 
sity  both  at  the  nadir  and  the  limb  to  profile  parameters  not  only  in  the  layer  peak  but  also  in  the 
topside  and  bottomside.  These  coefficients  are  simple  in  nature,  allowing  a  rapid  converion  op¬ 
erationally.  An  error  analysis  was  included,  de«;ribing  the  anticipated  uncertainties  in  the  results, 
from  (1)  internal  scatter  in  the  calibrations,  (2)  propagated  errors  from  the  errors  in  the  (nightglow) 
measurement,  and  (3)  uncertainties  in  the  reaction  rates,  expressed  as  partial  derivatives. 

A  complete  description  of  this  work  is  contained  in  the  report 

THE  DEVELOPMENT  OF  CALIBRATIONS  FOR  IONOSPHERIC  PARAMETERS 
FROM  DMSP-BASED  NIGHTGLOW  OBSERVATIONS  by  M.W.  Fox.  D.N.  Anderson  and 
R.E.  Daniell  (1993) 

a  draft  copy  of  which  has  been  delivered  to  PL/GPIM. 

2.5  A  Near-Global  Model  of  Total  Content 
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There  is  still  a  large  body  of  the  ionospheric  user  community  that  relies  on  simple  empirically- 
derived  descriptions  of  the  ionosphere.  A  number  of  these  models  use  coefficient-based  descrip¬ 
tions  of  peak  parameters  to  scale  a  profile  of  a  particular  shape,  and  TEC  is  derived  in  turn  from 
the  profile.  For  this  effort,  a  new  model  for  estimating  TEC  has  been  developed,  along  the  lines 
of  that  proposed  by  Brown  et  al.  (1991).  The  model  is  denoted  LMTEC,  denoting  Low  and  Mid¬ 
dle  latitude  TEC.  The  basic  premise  is  that  the  variations  of  slab  thickness  are  less  than  those  of 
foF2,  and  that  much  of  the  complexity  of  estimating  TEC  arises  from  the  estimate  of  foF2.  Total 
content  can  thus  be  estimated  simply  from  the  coefficient-based  (or  observed)  values  of  foF2  and 
a  model-based  estimate  of  slab  thickness. 

Firstly,  a  slab  thickness  model  needs  to  be  established  at  lower  latitudes.  At  middle  latitudes,  the 
model  of  Fox  et  al.  (1991)  has  been  used.  Also,  the  variations  of  slab  thickness  at  lower  latitudes 
in  the  Anderson  ( 1973)  model  have  been  described  in  the  same  manner,  and  the  latitude  variations 
of  these  coefficients  in  turn  summarized  numerically.  This  slab  thickness  model  has  been  delivered 
to  PIVGPIM. 

The  next  step  undertaken  was  validation  of  the  model  against  both  hourly  and  monthly  mean  TEC 
from  a  number  of  sites,  with  the  results  compared  with  those  achieved  by  other  models.  The  predic¬ 
tions  for  various  models  of  TEC  are  compared  with  the  monthly  mean  TEC  observed  from  Ramey. 
PR,  in  Figure  7.  Generally,  comparable  results  were  achieved,  though  overall,  the  LMTEC  results 
were  systematically  the  most  reliable.  This  was  most  evident  when  uniform  peak  parameters  had 
been  used  to  drive  the  various  models.  It  was  apparent  that  a  significant  source  of  the  numerical 
discrepancies  comes  from  not  precisely  estimating  the  times  of  the  post-sunrise  increase  and  the  di¬ 
urnal  maximum.  Possibly  some  improvement  would  be  obtained  in  LMTEC  if  higher  order  terms 
were  introduced  into  the  coefficients  that  describe  the  variations  of  slab  thickness.  Secondly,  a 
limiting  factor  at  lower  latitudes  in  using  LMTEC  will  be  the  accuracy  of  the  slab  thickness  model 
derived  from  the  Anderson  ( 1973)  theoretical  ionospheric  model.  The  results  suggest  that  better  re¬ 
sults  may  be  achieved  by  using  higher  order  terms,  both  in  the  Local  Time  and  latitudinal  variations 
of  the  model  slab  thickness. 

This  work  is  to  be  presented  as 

A  NEAR-GLOBAL  MODEL  OF  TOTAL  ELECTRON  CONTENT,  by  M.W.  Fox  (1993) 

at  the  Ionospheric  Effects  Symposium,  May  4-6,  1993,  in  Alexandria,  VA,  and  will  subsequently 
be  submitted  to  Radio  Science. 

2.6  A  New  Simple  Electron  Density  Profile  Description 
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A  variety  of  the  electron  density  profile  models  currently  in  use  rely  on  specification  of  the  layer 
peak  from  some  external  means  (e.g.  coefficients)  and  some  parameterized  layer  shi^.  How¬ 
ever,  a  recent  comparison  of  a  number  of  ionospheric  models  currently  in  use  (Brown  et  al.,  1991) 
showed  that  no  model  provided  a  really  good  description  of  the  layer  shape,  as  measured  by  the 
slab  thickness. 

For  this  effort,  a  new  expression  for  electron  density  profile  shape  has  been  developed.  The  system 
scales  the  profile  not  only  to  the  peak  but  also  to  the  value  of  slab  thickness  estimated  from  the 
model  described  in  the  previous  section  (note  also  that  total  content,  or  a  topside  density  could 
equally  well  be  used).  This  additional  requirement  is  necessary  to  constrain  the  profile  shape.  This 
approach  is  similar  to  the  SLIM  model  of  Anderson  et  al.  (1987)  in  that  it  consists  of  Chapman 
layers  in  the  topside  and  bottomside.  However,  in  this  case,  the  scale  height  of  the  Chapman  layers 
are  allowed  to  vary  linearly,  and  this  gives  rise  to  a  relatively  simple  (analytic)  expression  far  the 
layer  shape.  The  name  of  this  model,  DVCHAP,  refers  to  Double  Variable  (-scale  l^ight)  Chapman 
layers. 

The  key  to  specifying  the  ionosphere  in  DVCHAP  is  in  knowing  the  variations  of  the  parameters 
that  describe  scale  height  variations  in  the  topside  and  bottomside.  Two  studies  were  undertaken  to 
examine  this.  Firstly,  the  SLIM  approach  was  adopted,  and  the  variations  of  these  parameters  in  the 
Anderson  (1973)  theoretical  ionospheric  model  were  described.  Secondly,  the  Millstone  IS  Radar 
database  of  Buonsanto  (1989)  was  examined,  and  profile  residuals  were  studied  as  a  function  of 
how  many  ionospheric  parameters  were  known  beforehand.  Overall,  the  model  is  quite  successful 
in  describing  both  theoretical  and  observation  profile  shapes.  Examples  of  fits  to  Millstone  Hill 
profiles  are  seen  in  Figure  8. 

This  work  is  currently  being  prepared  for  publication  as 

DVCHAP:  A  NEW  SIMPLE  ELECTRON  DENSITY  PROFILE  MODEL,  by  M.W.  Fox 
(1993) 

and  will  be  submitted  to  Radio  Science. 

2.7  SKYMAPS  :  A  System  For  Calibrating  7774A  Images  And  TEC 

For  this  effort,  we  examined  the  feasability  of  a  system  that  would  permit  both  the  routine  moni¬ 
toring  of  total  electron  content  from  sites  with  all-sky  camera  facilities,  and  that  would  ultimately 
provide  a  means  of  establishing  a  calibration  for  sky-maps  (local  line-of-sight  variations)  of  total 
content  from  the  ground-receiver  sites  in  the  TISS  networic.  All-sky  7774X  images  obtained  at 
Millstone  Hill,  MA,  using  the  NSF/CEDAR  all-sky  imaging  system  run  by  Boston  University  (see 
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e.g.  page  258  of  Liu,  1989)  have  been  used  as  a  basis. 

The  questions  to  answer  are 

1)  what  are  the  variations  that  occur  over  any  given  site? 

2)  how  can  the  TEC  variations  best  be  recognized  against  a  background  that  includes  both  terrestrial 
and  non-terrestrial  sources  of  7774X  radiation  (or  how  can  the  system  best  be  trained)? 

3)  what  are  the  best  model  summaries  of  the  observed  features  so  that  recognition  of  a  feature  can 
lead  to  an  update  of  the  “skymap”? 

Initially,  a  database  of  these  observations  covering  a  range  of  geomagnetic  conditions  was  estab¬ 
lished.  One  criterion  for  the  level  of  “disturbance”  in  given  image  was  seen  to  be  the  amount  of 
scatter  about  an  average  intensity  vs.  sec(x)  relationship.  It  was  next  established  that  in  model 
ionospheres  that  a  reliable  calibration  exists  between  S  Ne  (w  TEC)  and  D  (or  F-region  O'*' 
radiative  recombination).  Comparisons  then  revealed  that  the  Bent  model  (LLewellyn  and  Bent, 
1973)  produced  a  model  ionosphere  that  was  most  consistent  with  the  7774  X  observations.  This 
model  also  has  the  advantage  of  being  updated  by  available  data,  so  would  be  more  appropriate  to 
the  proposed  ITSS-based  studies.  Unfortunately,  in  the  database,  the  disturbances  seen  in  the  im¬ 
ages  were  auroral  in  nature,  and  no  trough  passages  could  be  detected.  For  this  reason,  “training” 
of  the  system  is  waiting  on  a  larger  database  when  F-region  effects  will  be  included. 


11 


REFERENCES 


Anderson.  D.  N..  A  thewetkal  study  of  the  ionospheric  F  region  equatorial  anomaly,  I.  Theory. 
Planet.  Space  ScL.  21, 409, 1973. 

Anderson,  D.N..  Mendillo.  M.  and  Hemiter,  B.  A  semi-empirical  low-latitude  ionos|rfieric  model. 
Radio  Science,  22, 292, 1987. 

Brown,  L.D.,  Daniell,  R.E..  Fox,  M.W.,  Klobuchar,  J.A.  and  Doherty,  P.H.  Evaluation  of  six  iono¬ 
spheric  models  as  predictors  of  total  electron  content.  Radio  Science,  26,  p.  1007, 1991. 

Buonsanto.  M.J..  Comparison  of  incoherent  scatter  observations  of  electron  density,  and  electron 
and  ion  temperature  at  Millstone  Hill  with  the  International  Reference  Ionosphere,  J.  Atmos. 
Terr.  Phys.,  51, 441. 1989. 

Crowley.  G..  Emery.  B.A..  Roble,  R.G.,  Carlson,  H.C..  Salah,  J.E..  Wickwar,  V.B..  Miller,  K.L., 
Oliver,  W.L..  Burnside,  R.G.,  and  F.A.  Marcos,  Thermospheric  dynamics  during  September 
18-19, 1984;  2.  VaHdation  of  the  NCAR-TGCM,  J.  Geophys.  Res.,  94, 16945-16959, 1989. 

Daniell,  R.E.,  W.G.  Whartenby,  L.D.  Brown,  D.N.  Anderson.  M.W.  Fox,  P.H.  Doherty,  D.T.  Decker, 
J.J.  Sojka  and  R.W.  Schunk,  PRISM:  A  Parameterized  Real-time  Ionospheric  Specification 
Model,  to  be  submitted  to  Radio  Science,  1993. 

Fox,  M.W.,  The  Day-to-Day  Variability  of  Total  Content,  Peak  Density  and  Slab  Thickness,  and 
the  Ionospheric  Response  to  Geomagnetic  Storms,  GL-TR-90-0313,  Geophysics  Laboratory, 
Hanscom  AFB,  Bedford,  MA,  U.S.A.,  1990.  ADA233294 

Fox,  M.W.,  Mendillo,  M.  and  Klobuchar,  J.A.  Ionospheric  equivalent  slab  thickness  and  its  mod¬ 
eling  applications.  Radio  Science,  26,  p.  429, 1991. 

Fox,  M.W.,  D.N.  Anderson  and  R.E.  Daniell,  Developing  a  system  f  calibrations  for  ionospheric 
parameters  from  DMSP-based  nightglow  observations,  to  be  submitted  to  Radio  Science, 
1993. 

Hedin,  A.  E.,  MSIS-86  thermospheric  model,  J.  Geophys.  Res.,  92, 4649, 1987. 

Hedin,  A.  E..  Extension  of  the  MSIS  thermosphere  model  into  the  middle  and  lower  atmosphere, 
J.  Geophys.  Res.,  96, 1159-1172, 1991. 

Llewellyn,  S.K.  and  Bent,  R.B.  Documentation  and  Description  of  the  Bent  Ionospheric  Model, 


13 


APCRL-TR-73-06S7,  Air  Foice  Cambridge  Research  Lab..  Bedford.  MA.  U.S.A..  1973. 
AD772733. 

Liu.  C.H.  (Ed.)  World  lonosphere/Thermosphere  Study  (WITS)  HandbotA.  Volume  2.  Experimen¬ 
tal  Techniques.  December  1989. 

Pi.  X.Q..  M.  Mendillo.  M.W.  Fox  and  D.N.  Anderson.  Diurnal  double  maxima  patterns  in  the  F  re¬ 
gion  ionosphere:  Substorm-related  aspects.  Accepted  for  publication.  /  Geophys.  Res,  1993. 

Rawer,  K.  and  Piggott.  W.R.  (Eds.)  The  development  of  IRI-90  (Proceedings  of  woiksht^  held  at 
Abingdon,  U.K.),  Adv.  Space  Res.,  10  (11),  1989. 

Schoendorf,  J.,  Crowley.  G.,  and  R.G.  Roble.  Neutral  density  structures  in  the  high  latitude  ther¬ 
mosphere,  EOS  Trans.,  72, 353,  1991. 


H 


ALTITUDE  (KHI 


%TIOAL  DENSITY  PERTURBATION  (LAT.  0) 


II 

o 


ru 


H  <o 


H  00 


(Hxiaaniiiiv 

Figure  2 


16 


LOCAL  TINE  IHRSI 


LOCI*  ocNSirr  icn/cnsi  LO(ii«  ocmsiit  icm/chsi 


DAYTIME  DENSITY 


(mV/m) 


o  o 
•  • 

04  »- 


O  O 
•  • 

04 

I  I 


G"/!*  ,,0l)  031 


(o»s/ui)  'a 


20 


Fi^iire 


1001 


I 

I 


Romey  (M«on)  :  Ot>iefv«d  ond  Predicted  TEC  :  D«(oult  1oF2 
I 
I 


I 


Ramey  (Meon)  :  Observed  and  Predicted  TEC  ;  Default  (or2 


September  1983  October  1983  November  1983  December  1983 


Figure  7 


21 


Altitude  Altitude  Altitude 


10  September  1980  10  September  1980 


20  December  1980  20  December  1980 


17  Jonuory  1980  16  Jonuory  1980 


Figure  6 
7? 


APPENDIX 


FORTRAN  Codes  for  GMSIS83.  DCORR,  and  TMOD  Programs 


23 


o  o 


soBRoorm  nioD 


•ubroutln*  taK>d(apflO,itlaM,poa,ss,d«is,tM4>,pr«SB) 

c  calculata  grovas/msia  danaitlas  with  "corractad"  tidal  parturbations 

c  batwaan  80  and  145  kmi  aBooth  BMrgar  with  grovaa/naia  batwaan 

c  145  and  165  loa. 

raal*8  dana(7,50),  taap(50),  praaa(50) 
raal*4  apfl0(3),  zz(3),  poa(3) 
intagar*4  itina(5) 
dinanaion  am(2),  ph(2) 

c  Thia  aubroutina  ia  only  antarad  if  z.gt.80  and  z. It. 165 

c  It  ia  aaaunad  that  zz( l)*zz(2)"altituda>z 

xl«poa(l} 

t<Bitiiiia(4) 

iiio»itima(2) 

z-zz(l) 

if (z. la. 145. )  go  to  99 

c  If  z.gt.l45  intarpolata  with  grovaa/maia  at  165  km 

c  1.  Calculata  tidal  corraction  to  grovaa/maia  at  145  km 

call  dcorr(mo,xl,145.,t,am,ph,rcl) 

c  2.  Calculata  grovaa/maia  tidal  corraction  at  165  km 

zz(l)-165. 
zz(2)*165. 
call  aatawch ( 7 , 0 . ) 
call  aatawch (8,0.) 

call  grovaa ( apf 10 , itime, poa , zz , dana , tan^, preaa ) 
dmaan*dana (1,1) 
call  aatawch(7,l. ) 
call  aatawch ( 8, 1. ) 

call  grovaa (apflO,itlma, poa, zz, dana, tang), praaa) 

dt*dana(l,l) 

rc2-dt/dnMan 

c  3.  Linaarly  intarpolata  to  gat  tidal  danaity  correction  at  altitude  z 

rc-rcl+(rc2-rcl)*(z-145. )/20. 
go  to  100 

99  continue 

call  dcorr(mo,xl,z,t,am,ph,rc) 

100  call  aata%rch(7,0. ) 
call  aata%rch  (8,0.) 

call  grovaa ( apf 10 , itima , poa , zz , dana , tamp, praaa ) 
dana  ( 1 , 1  )«>dana  (1,1)  *rc 

return 

and 
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o  o 


PROGRAM  GMSIS83 


program  gmals83 

paraaMtar(nt*13,  ni«2S) 

raal*8  dana(7,50),t«i9(50),prass(50) 

raal*4  apfl0(3),zz(3),poa(3) 

lnt«gar*4  itl]M(S) 

dijaanslon  am(2),  ph(2) 

dlnansion  dt(nt) ,dr(nt,nz) ,it(nt) 

apfl0(l)-120. 

apfl0(2)-120. 

apfl0(3)-4. 

ltijBa(l)«>80 

ltiBM(2)*3 

itina(3)-21 

ltijM(4)"0 

itliM(5)-0 

zz(3)-1.0 

pos(2)«0.0 

c***  sat  swltchas  to  turn  off  semiannual,  annual,  UT,  J^,  longitude,  and 
c  terdlurnal  variations,  for  comparison  with  annual  mean  density 
c  perturbations  from  DCORR. 

call  sets%ich(3,0. ) 
call  setswch(4,0. ) 
call  sets%fch  (5,0.) 
call  sets%fch(6,0. ) 
call  setswch(9,0. ) 
call  setswch(10,0. ) 
call  sets%rch(ll,0. ) 
call  setswch (12,0.) 
call  sets%rch(13,0. ) 
call  setswch(14,0. ) 
c 

c  Switches  7  and  8  are  turned  on/off  later  to  turn  on/off  diurnal  and 
c  semidiurnal  tidal  variations  In  groves/msls. 
c 

c  ***  The  following  code  calculates  the  percent  total  mass  density  variation 
c  about  the  diurnal  mean  values  at  various  heights  and  latitudes  using  the 
c  groves/msls  model.  Tables  of  percent  deviation  vs.  height  (80-200  Icm)  and 
c  local  time  (0-24  hrs)  are  generated  at  latitudes  0,  30,  and  60  degrees 
c  latitude.  Other  Input  parameters  for  the  groves  subroutine  are  given  above, 
c 

c  Further  below,  a  similar  code  appears,  except  that  the  percent  total  mass 
c  density  variation  Is  provided  by  the  DCORR  subroutine  developed  by  Forbes, 
c  This  latter  code  Illustrates  the  liig>lementatlon  of  DCCXtR,  ^Ich  Is 
c  only  valid  between  80  and  150  km  (hence  the  tidal  switches  In  groves/msls 
c  must  be  turned  on/off  depending  whether  one  Is  outslde/lnslde  the  80-150 
c  km  height  regime.  In  the  Isplementatlon  shown  here,  DCORR  Is  actually  called 
c  within  subroutine  TMOD.  For  altitudes  between  145  and  165  km,  TMOD  provides 
c  a  linear  Interpolation  between  the  DCORR  results  at  145  km,  and  the 
c  groves/msls  results  at  165  km,  so  that  no  density  discontinuity  exists  as  a 
c  result  of  using  DCORR. 
c 

c  It  Is  Important  to  note  that  DCORR  only  corrects  the 

c  tidal  DEHSITY  variation;  that  Is  ,  when  It  Is  Implemented  In  the  following 
c  fashion,  as  a  correction  to  the  diurnal  mean  density  from  groves/msls,  that 
c  other  parameters  (tesq),  press,  etc.)  given  by  groves/msls  are  the  diurnal 
c  mean  values  with  NO  tidal  variation  whatsoever.  If  such  variations  are 
c  required,  the  groves/msls  model  would  have  to  be  called  again  with  the 
c  appropriate  Input  parameters. 
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e - - — — — — — — — - - - GROVBt  V«rsioB  —————— — — — 

300  fonMt(20x,'PBRCBIIT  TZML  OBHSITT  VARZATIOM  FROM  GIIOVBS/N8I803 '  ) 


print  300 


do  1-1,3 

xl-(l-l)*30. 

poa(l)-xl 

100  format ( lx, Ohaonth-, 13, 3x,4hlat-, f4 .0 ) 

print  100,  ltiaM(2),xl 
200  fonMt(lx,3hALT) 

print  200 

do  1-1, ns 
2-80.4(1-1) *5. 
sz(l)-s 
2S(2)-t 

davg-0 . 0 
do  j-l,nt 
itiM(4)-(  j-l)*2. 
it(j)-itiiaa(4) 
call  satswch(7,l. ) 
call  sats%rch(8,l. ) 

call  grovoa  ( apf  1 0 ,  itlna ,  poa ,  sz ,  dana ,  tanp ,  praaa ) 

dr(j,i)-dana(l,l) 

davg-davg4dana (1,1) 

and  do 

davg-davg-dana (1,1) 
davg-davg/12. 
do  j-l,nt 
dd-dr( j,l) 

j  *  1  )*100 .  *  (dd-davg)  /davg 
and  do 
and  do 

call  print(dr,nt,nz,lt) 
and  do 

c - 

c  ****  Tha  "aaaa”  coda  ia  givan  balow,  axcapt  DCORR  ia  uaad  to  provida  tha 
c  tidal  danaity  variationa.  Llnaa  of  coda  diffarant  froa  tha  abova  coda  ara 
c  givan  in  UPPERCASE  LETTERS. 


c - 

c - DCORR  Varaion - 

400  format (2 Ox, 'PERCENT  TIDAL  DENSITY  VARIATION  FROM  DCORR') 


print  400 

do  1-1,3 

xl-(l-l)*30. 

poa(l)-xl 

print  100,  itlma(2),xl 
print  200 

do  i-l,nz 
z-80.4(i-l)*5. 

davg-0 . 0 
do  j-l,nt 
iti]M(4)-(  j>l)*2. 
it(j)-itima(4) 
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»*(!)-» 

ss(2)-s 

ZF(I.OT.80.AHD.I.LT.165.)  TBBH 

CALL  TIIOD(APF10,ITZNB,PO8, 18, DBM8,TBMP, PRESS) 

ELSE 

CALL  8ET8IPCB(7,1.) 

CALL  SBTSliCBCS,!.) 

CALL  GROVE8(APF10,ITIME,PO8,lS, DEES, TEMP, PRESS) 
ENDIF 

dr(j,i)-d*n»(l,l) 

davg-davg4^ns  (1,1) 

•nd  do 

davg^avg-><ten«  (1,1) 
dav9*davg/12 . 
do  j-l,nt 
dd-dr( j,l) 

dr ( j , 1 )>100 . * (dd-davg ) /davg 
•nd  do 
•nd  do 

call  prlnt(dr,nt,nz,lt) 

•nd  do 


c - - - 

•top 

•nd 

•ubroutinn  prlnt(dr,nt,nz,it) 
dlJMnslon  dr(nt,nz),lt(nt) 

100  fonMit<lx,f4.0,2x,13f5.1) 

200  £or8Utt(lx,4hLST->,2x,  13(lx,12,2x) ) 

300  £oz«at(lBl) 

do  l>il,nz 
z-80.-f(l-l)*5. 

print  100,z,(dr( j,i), j«l,nt) 

•nd  do 

print  200,  (it(i),i«l,nt) 
print  300 

rotum 

•nd 
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e  suBRooTim  ocom 

c  _ 


subroutln*  dcorr(BO,xl,s,t,«fli,ph,rc) 

dlMiislon  ar(2,7),«l(2,7),p(7),xr(2),xi(2),aB(2),ph(2) 

conplM  dll,d22,d23,d24,d2S 

coav>l«x  fs(4},fd(4),B(4,12)/d(12},xd,XB 

coBBK>n /co«  f  f  B /xr ,  xl 

c  Coaq;>ut«B  dmBlty  corraction  factor  danx  -  danB*rc  to  account 

c  for  danaity  parturbationx  diM  to  tldaa  in  tha  SO  to  150  km 

c  halght  raglon.  Corractlon  factor  la  datarainad  froa  Lagandra 

c  Polynoodal  raconatructlon  of  diurnal  (n»l)  and  Baaidiumal  (n->2) 

c  harmonica  of  perturbation  danaity.  Thia  aubroutina  alao  raturna 

c  ralativa  danaity  anplituda  (daltarho/rho)  and  phaaa  (local  tlma 

c  of  maximum)  for  tha  n-1  and  n«2  harmonica. 

c  Thia  J^rll  1992  varalon  includaa  obBarvationally<Klatarminad  calibrating 

c  coafficlanta  for  tha  HMBa,  for  month-to-B»nth  variationa  in  tha 

c  n<*l  and  n*2  haraxinica. 

c  mo  -  month 

c  xl  -  latitude  in  dagraaa 

c  z  halght  in  km 

c  t  ■  local  aolar  tixM 

c  am  •  amplitude  of  ralativa  danaity  perturbation  at  (mo,xl,z) 

c  ph  "  local  tlma  of  maximum  of  ralativa  danaity  perturbation 

c  at  (ao,xl,s) 

c  n  ■  1  ia  tha  diurnal  harmonic,  with  24-hour  period 

c  n  ■  2  ia  tha  aamidiumal  hanaonic,  with  12-hour  period 

c  rc  •  danaity  corractlon  factor  at  (xl,z,t)  which  can  be  applied  to 

c  any  u|k>«^  atmoaphara  danaity  andal  which  doaa  not  include  tidal 

c  af facta  between  80  and  150  km.  In  acme  caaaa,  if  a  tidal 

c  affect  ia  included  in  auch  a  modal  above  aoma  height  zo,  than 

c  a  Bchama  muat  be  introduced  to  provide  a  amooth  tranaition  between 

c  tha  two  apaclficatlona  of  danaity  variation.  Thia  ahould  be  dona 

c  external  to  tha  praaant  atibroutina. 


c***  hough  function  valuaa  at  normalizing  latitudaa  0,  24,  36,  42 
c  dagraaa  for  (2,2),  (2,3),  (2,4),  (2,5)  modaa,  raapactivaly  *** 
data  h22,h23,h24,h25/l. 165, 1.091, 1.056, 1.055/ 
data  hll/1.0/ 


c  **  HOTB  ** 

C  THE  POLLOWING  TABLES  OF  COEFFICIENTS  ARE  THE  DRIVERS  FOR  THE  WHOLE  DENSITY 
c  CORRECTION  SCHEME  FOR  THE  80  -  150  XN  REGION.  THESE  TABLES  CAN  AND  SHOULD 
C  BE  UPDATED  AS  N<»E  DATA  BECOMES  AVAILABLE.  OPERATIONALLY,  THE  TIDAL  DENSITY 
c  C(»RECTION  ALGORITH  SHOULD  BE  IMPLEMENTED  SO  THAT  THE  FOLLOWING  COEEFICIBNTS 
C  CAN  BE  EASILY  UPDATED. 

c  ***  aamidlurnal  tmqpa  for  (2,2),  (2,3),  (2,4),  (2,5)  at  100  km 
c  from  Forbaa  and  Vial  (1989)  va.  month.  Thaaa  rapraaant  a  purely 

c  theoretical  aatimata  of  tha  coafficlanta  to  calibrate  tha  hmaa'a  *** 

c  data  ( (B(i,k), 1-1,4), k«l, 12)/ 

c  1  (2.76,0.5),  (1.83,9.3),  (3.77,1.4),  (7.03,11.4), 
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c 

2  (2.18,0.7), 

(2.17,8.7), 

(4.14,1.1), 

(6.17,11.6), 

e 

3  (1.27,10.8), 

(.416,6.5), 

(9.28,1.3), 

(.402,3.4), 

c 

4  (1.63,0.2), 

(.817,6.4), 

(8.11,1.2), 

(3.52,6.7), 

c 

5  (1.87,0.4), 

(.932,1.8), 

(1.35,0.9), 

(4.30,6.1), 

c 

6  (2.09,0.4), 

(1.49,1.7), 

(.996,3.1), 

(3.40,6.0), 

c 

7  (1.97,0.7), 

(1.85,0.1), 

(2.60,6.0), 

(4.70,5.8), 

c 

8  (2.10,0.7), 

(2.47,11.7), 

(2.12,11.7), 

(5.34,6.2), 

e 

9  (1.22,11.2), 

(.473,8.8), 

(4.52,0.6), 

(3.93,7.0), 

c 

1  (1.64,0.3), 

(.965,4.3), 

(6.66,1.3), 

(1.14,2.9), 

c 

2  (2.50,0.3), 

(1.29,7.7), 

(4.09,0.9), 

(5.24,11.9), 

c 

3  (2.66,0.3), 

(2.32,8.5), 

(3.03,1.8), 

(5.36,11.6)/ 

c  ***  •Mddlurnal  taapa  for  (2,2),  (2,3),  (2,4),  (2,5)  at  100  Im 
c  froai  fit  to  Monthly  cllMatologlcal  Nl/r  radar  data  80-110  ka; 
c  Incoharant  acattar  radar  data  not  Ineludad  in  flta  ********* 


c 


c 

data  ((s(i,k),l 

-l,4),k-l,12)/ 

c 

1  (1.97,1.3), 

(2.51,2.0), 

(1.58,9.3), 

(2.24,10.7), 

c 

2  (1.56,0.3), 

(1.14,3.3), 

(  .72,10.2), 

(2.12,11.5), 

c 

3  (1.55,1.2), 

(  .87,10.2), 

(1.83,11.8), 

(  .30,1.6), 

c 

4  (3.07,1.9), 

(2.19,10.0), 

(1.34,11.4), 

(2.15,5.6), 

c 

5  (1.31,1.8), 

(3.07,9.6), 

(1.55,  1.5), 

(1.56,5.4), 

c 

6  (1.58,2.0), 

(1.25,11.8), 

(1.02,10.8), 

(  .68,5.0), 

c 

7  (2.00,2.5), 

(1.37,10.5), 

(  .87 ,1.8), 

(  .23,7.0), 

c 

8  (2.53,2.9), 

(5.35,10.5), 

(1.18,10.5), 

(2.11,5.9), 

c 

9  (2.85,0.6), 

(4.28,9.4), 

(2.07,3.0), 

(2.00,6.3), 

c 

1  (2.09,2.0), 

(2.16,10.3), 

(  .90,10.2), 

( 1 . 06,5 . 7 ) , 

c  , 

2  (2.55,1.8), 

(  .80,3.2), 

(  .51,10.4), 

(  .54,10.4), 

c 

3  (2.29,2.3), 

(1.46,1.2), 

(2.46,10.8), 

(1.26,10.4)/ 

c  *** 

saaidiumal  tanps 

for  (2,2),  (2,3) 

,  (2,4),  (2,5) 

at  100  ka  froa 

c  fit  to  okonthly  eliaatologlcal  IS  *AIID*  MLT  radar  data  80-130  kn; 


c 


data  ((s(l,k) 

, i-1, 4 ),k-l, 12)/ 

1 

(1.29,5.3), 

(2.63,0.3), 

(2.00,10.7), 

(1.44,10.2), 

2 

(0.42,9.6), 

(1.01,10.9), 

(1.00,10.7), 

(1.08,11.8), 

3 

(1.00,10.4), 

(2.88,10.3), 

(2.39,11.5), 

(0.42,5.9), 

4 

(1.05,1.1), 

(4.58,10.3), 

(2.13,0.8), 

(2.61,5.3), 

5 

(1.12,0.1), 

(3.58,10.3), 

(1.52,1.2), 

(1.69,4.9), 

6 

(  .29,1.7), 

(2.76,10.6), 

(1.53,11.5), 

(1.62,4.9), 

7 

(  .61,3.1), 

(3.25,10.5), 

(1.78,1.1), 

(1.40,6.1), 

8 

(3.08,3.9), 

(5.18,10.5), 

(2.72,6.1), 

(5.18,10.5), 

9 

(1.50,11.0), 

(3.84,8.4), 

(3.52,11.2), 

(2.39,5.9), 

1 

(1.54,5.2), 

(3.63,9.5), 

(1.68,11.5), 

(1.93,5.0), 

2 

(1.34,4.3), 

(1.74,11.1), 

(1.49,0.0), 

(  .47,7.5), 

3 

(1.36,3.9), 

(2.28,11.3), 

(2.48,11.1), 

(  .45,10.0)/ 

c****  Diurnal  calibrating  coafficianta  for  (1,1)  Moda  datarmlnad  fr<Hi 
c  RaMy,  P.R.  and  Townasvilla,  Australia  wind  Maaauranants  at  90  ka. 
c  (anp,  i^asa)  Ai^lltuda  raprasants  tans  of  SMtars/sac  of  diurnal 

c  aastward  wind  aaplituda  at  90  ka  and  18  dagraas  latituda;  phasa  is 

c  local  tiaa  of  aaxiwua  eastward  wind.  k*l,12  corrasponds  to  jan,...dac. 

data  (d(k),k-l,12)/(2.0,12.),(2.5,12.),(4.0,13.),(2.5,13.), 

1  (2. 5, 13.), (1.5, 13.), (2. 0,14.), (4. 0,14.), (3. 0,14.), (1.0, 18.), 

2  (1.0, 8.), (2. 0,12.)/ 


C 


if (z. It. 80. or. 2. gt. 150.)  go  to  99 
ploacos(-l. ) 

c  Dafina  coaplax  HME's  at  (xl,s)t 
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emll  <llMMll(s,sl,dll) 
call  ?2(g,«l,d22) 
call  dh— ?2(t,»l,d23) 
call  dlMM24(s,sl,d24) 
call  dluMi25(s,xl,d25) 

c  Convart  noraalislng  factors  for  acmth  •  bk>  to  cosqplaxt 

do  1*1,4 
fs(i)*s(i,aK>) 
and  do 

call  fxbnd(fs,4,2) 

ld(l)-d(*o) 

call  fxbnd(fd,l,l) 

c  Dafina  diurnal  and  saaddlurnal  dansity  variation  (conqplax)  at 
C  (Z,Xl,IBO)t 

xd  -  hll*dll*fd(l) 

xs  -  -h22*d22*fs(l)  >  h23*d23*fs(2) 

6  -h24*d24*fs(3)  -  h25*d25*fs(4) 

xr(l)*raal(xd) 

xl(l)*alna9(xd) 

xr(2)*raal(xs) 

xl(2)*alna9(xs) 

c  Convart  from  coB9>lax  to  aB9>lltuda  and  phasa 
call  ap(xr,xl,aB,ph) 

c  Conputa  dansity  corractlon  factor  at  (xl,s,t) 

tc»l. 
do  n*l,2 

rc  -  rc  +  .01*aii(n)*cos(  (t>ph(n)  )*float(n)*pl/12. ) 
and  do 

99  ratum 

and 

subroutlna  ap(xr,xl,am,ph) 
dlissnolon  xr(2)/Xl(2),aiB(2),ph(2) 

c  convarts  coo^lax  fora  to  anplltuda  and  phasa  (hour  LST) 

pi  *  acos(-l.) 

<k>  n  *  1,2 

aa(n)  •  sqrt(xr(n}*xr(n)  xi(n)*xl(n)) 

If (aa(n) .aq.0.0)  than 

ih(n)  *  0.0 

alsa 

ph(n)  *  atan2(xl(n),xr(n) )*12./(pl*float(n) ) 
and  If 

If (ph(n) .lt.0.0)  ph(n)  *  ph(n)  *  24./float(n) 
and  do 

ratum 

and 


subroutlna  fxbnd(x,lc,n) 

c  convarts  aaplituda/phasa  (LT  of  aaxlaua)  to  coaplax  fora 
coaplax  x(4) 
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fn*float(n) 
pl*acos(-l . ) 


do  1  i-l,lc 
•B^r««l(x(l) ) 
phz">aJjuig(x(i) ) 
ang-phz*fn*pi/12 . 
xr-cos(ang) 
xi-sln(ang) 

1  x(i)*wip*cng>lx(xr,xl) 

roturn 

•nd 

•ubroutlno  dl«g(xl,p) 
dlownslon  p(7),dp(7) 
doublo  proclaion  dp,dxl,cs,fl 

c  conputos  ZiOgondro  polynmnlal  toma  through  ordar  6  at  latituda  -  xl 

cxl-90.-xl 
dxl  -  cxl 

ca  -  dcoa(.01745329252*dxl) 
dp(l)  -  1. 
dp(2)  •  ca 
p(l)  -  1. 
p(2)  -  dp(2) 
do  1-2,6 
£1-1 

dp(l+l)-((2.*£l-l.)*ca*dp(l)-(£l-l.)*dP(i-l))/fl 

p(l+l)  -  dp(l+l) 
and  do 

raturn 

and 


aubroutlna  dhBMll(z,xl,dll) 
cooqplax  dll 

dliaanalon  cr(15,7),cl(15,7),ar(7),al(7),p(7) 

C  Thin  aubroutlna  Intarpolataa  tablaa  of  6th-ordar  Ziagandra  polynomial 
C  coafflclanta  (plua  izaan)  (cr-raal;  cl-lnag),  z-  80-150  km  (5  km)  to 
C  glva  coafflclanta  ar,al  at  halght  -  z  (km)  and  latituda  -  xl  (dag) 

C  of  tha  con^lax  HME  for  parturbatlon  danalty  (parcant)  of  tha  dll  moda. 

data  ((cr(j,k),k-l,7),(cl(j,k),k-l,7),  j-1,5)/ 

6  0.055742,  0.008325,-0.361614,-0.006649,  0.659095,-0.006027,-0.434787 

6  ,-0.258680,  0.001143,  0.431332,-0.000285,-0.141635,-0.001914,-0.091767 

6  ,-0.058174,  0.007919,-0.349165,-0.006532,  1.004046,-0.005700,-0.763388 

6  ,-0.518292,  0.003793,  1.047003,-0.001672,-0.705835,-0.004603,  0.116936 

6  ,-0.539064,-0.001736,  1.172848,  0.001360,-0.931498,  0.001963,  0.275171 

6  ,-0.024488,  0.015533,  0.523313,-0.007622,-1.109352,-0.017836,  0.757128 

6  ,-0.080068,  0.013846,  0.533738,-0.009426,-0.883046,-0.013050,  0.501360 

6  ,  0.442837,  0.006633,-1.058604,-0.002871,  0.981687,-0.007695,-0.369514 

6  ,  0.253618,-0.012014,-0.425990,  0.007412,  0.157004,  0.011411,  0.062435 

6  ,  0.245779,  0.007735,-0.946197,-0.007577,  1.272626,-0.002400,-0.643730 


data  ((cr(j,k),k-l,7),(ci(j,k),k-l,7),  j-6,10)/ 

0.336366,-0.002470,-1.074346,-0.000243,  1.279186,  0.006014,-0.591427 
,-0.058599,  0.014701,-0.238977,-0.007744,  0.657508,-0.014746,-0.437246 
,  0.243605,  0.007443,-1.290724,-0.003386,  2.083051,-0.008559,-1.226405 
,-0.290373,  0.010932,  0.621476,-0.007397,-0.474814,-0.009338,  0.123319 
,-0.001579,  0.013160,-0.485392,-0.011644,  1.079454,-0.007430,-0.735050 
,-0.355180,  0.001790,  1.497327,  0.001020,-2.240973,-0.004378,  1.302434 
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,-0.137252,-0.012308,  0.488358,  0.011197,-0.804782,  0.006405,  0.308251 
,-0.227223,  0.008073,  1.387512,-0.005586,-2.403057,-0.003632,  1.502747 
,-0.173369,-0.003843,  0.983320,  0.000990,-1.644342,  0.003499,  1.000086 
,-0.061628,  0.006710,  0.735354,-0.002597,-1.453409,-0.007241,  0.959296 
/ 


data  ((cr(j,lc),k-l,7),(ci(j,k),lc-l,7),  j-11,15)/ 

-0.113343,  0.000173,  0.916851,  0.000168,-1.665328,-0.001011,  1.043190 
,  0.022661,  0.008209,  0.186905,-0.004531,-0.498586,-0.008897,  0.368341 
,-0.067469,  0.003385,  0.745517,-0.002250,-1.427778,-0.003014,  0.909452 
,  0.055857,  0.004795,-0.164810,-0.001071,  0.172236,-0.007848,-0.058354 
,-0.021591,  0.004696,  0.547862,-0.002947,-1.129912,-0.002504,  0.737553 
,  0.060800,  0.004317,-0.337485,-0.001998,  0.530214,-0.003967,-0.293422 
,  0.019602,  0.004867,  0.329300,-0.002868,-0.766190,-0.002282,  0.512607 
,  0.046687,  0.003427,-0.382124,-0.002847,  0.661418,-0.000227,-0.384560 
,  0.050399,  0.004739,  0.119047,-0.002744,-0.385148,-0.004707,  0.262498 
,  0.026298,-0.000052,-0.359160,-0.000394,  0.662306,  0.000642,-0.387582 
/ 


jz-((Z-80.)/5.  *  .001) 
dz«(z-80.-float( jz)*5. )/5. 

c  Modify  asymptotic  behavior  of  Legendre  polynomial  fit  at  high 
c  latitudes  to  eliminate  spurious  wiggles, 
f act*l . 0 
xm*abs(xl) 

if(xl.le.-70.)  xl— 70. 
if (xl.ge.-t-70. )  xl«+70. 

If (abs(xl) .ge.70. )  fact-exp(-.01*( (xl-xm)**2) ) 

call  dleg(xl,p) 

xr-0.0 

xl->0 . 0 

do  1  k«l,7 

ar(k)-cr( jz+l,k)*(l.-dz)  +  dz*cr( jz+2,k) 
ai(k)-ci(jz+l,k)*(l.-dz)  +  dz*ci( jz+2,k) 
xr-xr+ar ( k ) *p( k ) 
xl^xl-t-al  ( k )  *p(  k) 

1  continue 

dll*cmplx(xr,xi)*fact 

return 

end 


subroutine  dhme22(z,xl,d22) 
complex  d22 

dl^nslon  cr(15,7)  ,ci(15,7),ar(7),ai(7),p(7) 

C  This  subroutine  Interpolates  tables  of  6th-order  Legendre  polynomial 
C  coefficients  (plus  mean)  (cr-real;  cl-imag),  z-  80-150  km  (5  km)  to 
C  give  coefficients  ar,al  at  height  «  z  (km)  and  latitude  -  xl  (deg) 

C  of  the  con^lex  HME  for  perturbation  density  (percent)  of  the  d22  mode. 


data  ((cr(j 

,  k) ,k*l, 7) 

,(ci( j,k). 

k-1,7),  j-1,5)/ 

-0.169341, 

0.000000, 

0.264998, 

0.000000,-0.102197, 

0.000000, 

0.000000 

,-0.098691, 

0.000000, 

0.149398, 

0.000000,-0.053681, 

0.000000, 

0.001919 

,-0.361645, 

0.000000, 

0.475024, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,-0.206671, 

0.000000, 

0.315332, 

0.000000,-0.116397, 

0.000000, 

0.004485 

,-0.452214, 

0.000000, 

0.710608, 

0.000000,-0.286755, 

0.000000, 

0.000000 

,-0.235008, 

0.000000, 

0.349767, 

0.000000,-0.135104, 

0.000000, 

0.022193 

,-0.580654, 

0.000000, 

0.733212, 

0.000000,  0.000000, 

0.000000, 

0.000000 
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,-0.285256, 

0.000000, 

0.356361, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

,-0.683079, 

0.000000, 

1.049176, 

0.000000, 

-0.381796, 

0.000000, 

0.017647 

,-0.279810, 

/ 

0.000000, 

0.484218, 

0.000000, 

-0.288153, 

0.000000, 

0.085742 

data  ((cr(j 

,  k) , k— 1 , 7 ) 

.(ci< 

lt-1,7),  j- 

6,10)/ 

-0.785786, 

0.000000, 

1.053633, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

,-0.271038, 

0.000000, 

0.441329, 

0.000000, 

-0.181850, 

0.000000, 

0.003095 

,-0.916167, 

0.000000, 

1.176281, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

,-0.384130, 

0.000000, 

0.474877, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

,-1.180274, 

0.000001, 

1.461481, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

,-0.336436, 

0.000000, 

0.376645, 

0.000000, 

0.011233, 

0.000000, 

0.000000 

,-1.357902, 

-0.000055, 

1.637778, 

-0.000006, 

0.000000, 

0.000000, 

0.000000 

,-0.123959, 

-0.000028, 

0.066249, 

-0.000003, 

0.000000, 

0.000000, 

0.000000 

,-1.393752, 

0.000009, 

1.663134, 

0.000001, 

0.000000, 

0.000000, 

0.000000 

,  0.227522, 
/ 

0.000005, 

-0.376663, 

0.000001, 

0.000000, 

0.000000, 

0.000000 

data  ((cr(j 

,  k) ,k— 1,7) 

» (cl( j,k) , 

k-1,7),  j- 

11,15)/ 

-1.289872, 

-0.000001, 

1.575239, 

-0.000001, 

-0.211185, 

0.000003, 

-0.086052 

,  0.409935, 

-0.000001, 

-0.737090, 

-0.000001, 

0.301875, 

0.000001, 

0.034503 

,-1.162509, 

0.000000, 

1.362541, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

,  0.539921, 

0.000000, 

-0.915434, 

0.000000, 

0.391622, 

0.000000, 

0.002363 

,-1.008992, 

0.000000, 

1.170860, 

0.000000, 

-0.098578, 

0.000000, 

-0.075953 

,  0.669268, 

0.000000, 

-1.066391, 

0.000000, 

0.411981, 

0.000000, 

0.000000 

,-0.849756, 

0.000000, 

0.901213, 

0.000000, 

-0.078384, 

o.ocoooo. 

-0.055649 

,  0.777714, 

0.000000, 

-1.162062, 

0.000000, 

0.386092, 

0.000000, 

0.000000 

,-0.705545, 

0.000000, 

0.685605, 

0.000000, 

0.044003, 

0.000000, 

0.000000 

,  0.852396, 

0.000000, 

-1.061791, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

6  / 


j*-((*-80.)/5.  +  .001) 
d*-(*-80."£loat( j*)*5. )/5. 

c  Modify  asymptotic  behavior  of  Legendre  polynomial  fit  at  high 
c  latitudes  to  eliminate  spurious  wiggles. 
fact-1.0 
xm-abs(xl) 

If (xl. le.-70. )  X1--70. 

If  (xl.ge.'f70. )  xl— »-70. 

If (abs(xl) .ge.70. )  fact-expC^.Ol^C (xl-xm)**2) ) 

call  dleg(xl,p) 

xr-0.0 

xl-0 . 0 

do  1  k-1,7 

ar(k)-cr( jz<>-l,k)*( l.-dz)  dz*cr( jz-»'2,k) 
al(k)-cl(  jz-»-l,k)*(  l.-dz)  +  dz*ci(  jz-f2,k) 
xr-xr+ar ( k ) *p ( k ) 
xl-xi-i-ai(k)*p(k) 

1  continue 

d2 2-cmplx ( xr , xl ) * f act 

return 

end 


subroutine  dhme23( z,xl,d23) 
conqplex  d23 

dimension  cr(15,7) ,cl( 15,7),ar(7) ,ai(7),p(7) 

C  This  subroutine  Interpolates  tables  of  6th-order  Legendre  polynomial 
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C  coefficients  (plus  nean)  (cr«real|  ci^iMsg),  s*  80-150  kn  (5  ke)  to 
C  give  coefficients  ar,ai  at  height  -  s  (ka)  and  latitude  -•  xl  (d^) 

C  of  the  coaplex  HMB  for  perturbation  density  (percent)  of  the  d23  node. 


data  ((cr( j,k),k-l,7) 

,(ci( j,k). 

k-1,7),  j- 

1,5)/ 

0.000000,  0.006140, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

i 

0.000000,-0.011506, 

0.000000, 

0.034069, 

0.000000, 

-0.026685, 

0.000000 

t 

0.000000,  0.019873, 

0.000000, 

-0.027559, 

0.000000, 

0.010279, 

0.000000 

9 

0.000000,-0.015649, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.003407, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.073662, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,  0.022647, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.088268, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.189464, 

0.000000, 

0.183542, 

0.000000, 

0.000000, 

0.000000 

9 

/ 

0.000000,  0.032785, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

data  ((cr(j,k),k-l,7) 

1  (®i(  j  » 

k-1,7),  j- 

6,10)/ 

0.000000,-0.105329, 

0.000000, 

-0.087378, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,  0.013454, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.347701, 

0.000000, 

0.650321, 

0. 000000, • 

-0.410056, 

0.000000 

9 

0.000000,-0.432722, 

0.000000, 

0.869166, 

0.000000, 

-0.368488, 

0.000000 

9 

0.000000,-0.578832, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.195348, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.763669, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,  0.349028, 

0.000000, 

-0.785404, 

0.000000, 

0.477196, 

0.000000 

9 

0.000001,-0.711042, 

-0.000004, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

/ 

0.000004,  0.847518, 

-0.000007, 

-1.179172, 

0.000010, 

0.000000, 

0.000000 

data  ((cr(j,k),k-l,7) 

/(ci( 

k-1,7),  j- 

11,15)/ 

0.000024,-0.837956, 

-0.000058, 

0.974664, 

0.000084, 

0.007335, 

-0.000070 

9 

0.000031,  0.541760, 

-0.000082, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000000,-0.575385, 

0.000000, 

0.513577, 

0.000000, 

0.142813, 

0.000000 

9 

0.000000,  1.085389, 

0.000000, 

-1.418761, 

0.000001, 

0.000000, 

0.000000 

9 

0.000019,-0.334993, 

-0.000047, 

0.145217, 

0.000069, 

0.193542, 

-0.000057 

9 

0.000025,  0.665488, 

-0.000067, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000028,-0.143261, 

-0.000068, 

-0.075871, 

0.000100, 

0.155025, 

-0.000083 

9 

0.000035,  0.664996, 

-0.000094, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

9 

0.000026,-0.033387, 

-0.000064, 

-0.105135, 

0.000093, 

0.027952, 

-0.000077 

9 

0.000030,  0.598951, 

-0.000080, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

&  / 


jz-((r-80.)/S.  *  .001) 
d*-(z-80.-float( j2)*5. )/5. 

c  Modify  asymptotic  behavior  of  Legendre  polynomial  fit  at  high 
c  latitudes  to  eliminate  spurious  wiggles. 
fact«1.0 
xm*abs(xl) 

if (xl.le.-70.)  xl— 70. 
if  (xl.ge.-f70. )  xl«-f70. 

if (abs(xl} .ge.70. )  fact«exp(-.01*( (xl-xm)**2) ) 

call  dleg(xl,p) 
xr-0 . 0 
xl*0.0 
do  1  k-1,7 

ar(k)-cr( jz-fl,k)*( l.-dz)  -f  dz*cr( jz-f2,k) 
ai(k)>ci( jz-fl,k)*( l.-dz)  -f  dz*ci( jz-f2,k) 
xr-xr-f  ar  ( k )  *p(  k ) 
xl-xl-fal(k)*p(k) 

1  continue 

d2 3«cmplx ( xr ^ xi ) * f ac t 
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return 

•nd 


subroutine  dhne24(z,xl,d24) 
con^lex  d24 

dimension  cr(15,7),cl(lS,7),ar(7),al(7),p(7) 

C  This  subroutine  Interpolates  tables  of  6th-order  Legendre  polynomial 
C  coefficients  (plus  mean)  (cr-real|  clolmag),  z*  80-150  km  (5  km)  to 
C  give  coefficients  ar,al  at  height  -  s  (km)  and  latitude  -  xl  (deg) 

C  of  the  complex  HMB  for  perturbation  density  (percent)  of  the  d24  mode. 


data  ((cr(j 

,k),k-l,7) 

,(cl( j,k). 

k-1,7),  j-1,5)/ 

-0.001112, 

0.000000, 

-0.000112, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.003375, 

0.000000, 

0.008271, 

0.000000,-0.015943, 

0.000000, 

0.000000 

,-0.003972, 

0.000000, 

0.000654, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.005666, 

0.000000, 

0.011820, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.011253, 

0.000000, 

-0.031639, 

0.000000,  0.001761, 

0.000000, 

0.000000 

,  0.007032, 

0.000000, 

0.026286, 

0.000000,-0.109777, 

0.000000, 

0.068054 

,  0.009666, 

0.000000, 

0.004729, 

0.000000,-0.038965, 

0.000000, 

0.000000 

,  0.032353, 

0.000000, 

-0.013593, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.078647, 

0.000000, 

0.007455, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.053529, 
/ 

0.000000, 

0.176844, 

0.000000,-0.414525, 

0.000000, 

0.187836 

data  ((cr(j 

,  k) ,k*l, 7 ) 

,  (ci( j,k). 

k-1,7),  j-6,10)/ 

0.139092, 

0.000000, 

0.352783, 

0.000001  -0.936439, 

0.000001, 

0.469183 

,-0.085842, 

-0.000004, 

0.088246, 

-0.000005,  0.000000, 

0.000000, 

0.000000 

,-0.039680, 

-0.000012, 

0.035709, 

-0.000017,  0.000000, 

0.000000, 

0.000000 

,-0.188055, 

0.000028, 

-0.346149, 

0.000078,  0.859258, 

0.000073, 

0.000000 

,-0.210154, 

0.000000, 

-0.615333, 

0.000000,  1.138374, 

0.000000, 

0.000000 

,-0.068165, 

0.000000, 

0.018383, 

0.000001,-0.083173, 

0.000001, 

0.144534 

,-0.178775, 

-0.000001, 

-0.506718, 

-0.000001,  0.778054, 

0.000002, 

0.000000 

,  0.066324, 

0.000004, 

0.232282, 

0.000002,  0.000000, 

0.000000, 

0.000000 

,-0.142500, 

-0.000039, 

-0.113290, 

-0.000030,  0.149399, 

0.000069, 

0.094341 

,  0.103712, 
/ 

0.000056, 

0.379078, 

0.000006,  0.000000, 

0.000000, 

0.000000 

data  ((cr(j 

,k) ,k“l,7) 

,(ci( j,k). 

k-1,7),  j-ll,l5)/ 

-0.065776, 

0.000002, 

0.152067, 

0.000002,-0.428703,- 

-0.000003, 

0.347299 

,  0.095648, 

-0.000006, 

0.414288, 

-0.000001,  0.000000, 

0.000000, 

0.000000 

,-0.051485, 

0.000000, 

0.161603, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.103733, 

0.000000, 

0.335615, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,-0.031299, 

0.000000, 

0.202792, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.086790, 

0.000000, 

0.247737, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,-0.011905, 

0.000000, 

0.209203, 

0.000000,  0.000000, 

0.000000, 

0.000000 

,  0.086470, 

0.000000, 

0.411871, 

0.000000,-0.612001, 

0.000000,- 

-0.050880 

,  0.023170, 

0.000000, 

0.424157, 

0.000000,-0.739253, 

0.000000, 

0.198301 

,  0.077735, 

0.000000, 

0.324587, 

0.000000,-0.447950, 

0.000000, 

0.000000 

ft  / 


jz-((z-80.)/5.  -f  .001) 
dz*(z-80.-float( jz)*5. )/5. 

c  Modify  asymptotic  behavior  of  Legendre  polynomial  fit  at  high 
c  latitudes  to  eliminate  spurious  wiggles, 
f act>l . 0 
xm*abs(xl) 

if (xl.le.-70.)  X1--70. 
if (xl.ge.-t-70. )  xl>*'f70. 

if (abs(xl) .ge.70. )  fact»exp(-.012*( (xl-xm)**2) ) 


call  dlag(xl,p) 
xr»0.0 
xi«0.0 
do  1  k-1,7 

ar(k)*cr( jz'i'l,k)*(l.-^z)  ♦  dz*cr( jz’*-2,k) 
ai(k)-ci( jz-»-i,k)*(l.-dz)  -f  dz*ci(jz-*-2,k) 
xr-xr'i'ar(k)*p(k) 
xl-xl-t'al(k)*p(k) 

1  contlnua 

d2  4«cziplx  ( xr ,  xl )  *  f  act 

raturn 

and 


aubroutlna  dhiM25(z,xl,d2S) 
complax  d25 

dinanslon  cr(15,7),cl(15,7),ar(7),al(7),p(7) 

C  This  aubroutlna  Intarpolatas  tablas  of  6th-ordar  Lagandre  polynonlal 
C  coafficiants  (plus  izaan)  (cr«raal}  ci-lnag),  z«  80-150  km  (5  km)  to 
C  giva  coafficiants  ar,ai  at  haight  -  z  (km)  and  latituda  ■  xl  (dag) 

C  of  tha  conqplax  HMB  for  parturbation  dansity  (parcant)  of  tha  d25  mods. 


data  ((cr( j,k),k-l,7) 

,(ci( j,k). 

k-1,7),  j- 

1/5)/ 

0.000000,-0.026338, 

0.000000, 

-0.011444, 

0.000000, 

0.064545, 

0.000000 

i 

0.000000,-0.011902, 

0.000000, 

-0.007163, 

0.000000, 

0.000000, 

0.000000 

f 

0.000000,-0.055506, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,-0.027241, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,-0.123328, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,  0.059905, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,-0.039259, 

0.000000, 

-0.061414, 

0.000000, 

0.095220, 

0.000000 

0 

0.000000,  0.257230, 

0.000000, 

0.121375, 

0.000000, 

-0.549515, 

0.000000 

0 

0.000000,  0.339623, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

/ 

0.000000,  0.217038, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

data  ((cr( j,k),k-l,7) 

,(ci(j,k). 

Ic-l,7),  j- 

6,10)/ 

0.000000,  0.440520, 

0.000002, 

0.355003, 

0.000003, 

-0.882318, 

0.000007 

0' 

-0.000001,-0.176003, 

-0.000003, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

,-0.000024,-0.023348, 

-0.000067, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000003,-0.323842, 

0.000015, 

-0.311437, 

0.000029, 

0.753792, 

0.000064 

-0.000001,-0.195422, 

-0.000002, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,-0.069938, 

0.000001, 

0.106688, 

0.000002, 

0.000000, 

0.000000 

0 

0.000000,-0.115362, 

-0.000002, 

0.037199, 

-0.000003, 

0.000000, 

0.000000 

0 

0.000000,  0.057247, 

0.000000, 

0.602281, 

0.000001, 

-0.851582, 

0.000002 

0 

0.000000,-0.012605, 

0.000000, 

0.412166, 

0.000000, 

-0.593142, 

0.000000 

0 

1 

0.000000,  0.128748, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

data  ((cr( j,k),k-l,7) 

,(ci(j,k). 

Jc-l,7),  j- 

11,15)/ 

0.000000,  0.066114, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,  0.078923, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,  0.086642, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,  0.023933, 

0.000000, 

0.084672, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,-0.012868, 

0.000000, 

0.265220, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,  0.002520, 

0.000000, 

-0.005293, 

0.000000, 

0.125529, 

0.000000 

0 

0.000000,  0.065428, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,-0.001040, 

0.000000, 

0.019997, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,  0.079149, 

0.000000, 

0.000000, 

0.000000, 

0.000000, 

0.000000 

0 

0.000000,  0.006662, 

0.000000, 
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0.000000, 

0.000000, 

0.000000, 

0.000000 

&  / 


j*-({Z-80.)/5.  +  .001) 
d>-(s-80.-flo«t( js)*5. )/5. 

c  Modify  asyn^totic  bohavlor  of  Logondro  polynoodal  fit  at  high 
c  latltudas  to  aliainata  spurious  wigglas. 
fact-1.0 
xitt-abs(xl) 

if (xl.la.-70. )  X1--70. 
if  (xl.ga.<f70. )  xl*M‘70. 

if (abs(xl) .ga.70. )  fact-axp(-.01*( (xl-xm)**2) ) 

call  dlag(xl,p) 

xr-0.0 

xi-0 . 0 

do  1  k-1,7 

ar(k)-cr(  jZ'fl,k)*(l.-dz)  +  dz*cr(  jz+2,k) 
ai(k)-ci(  jz-»-l,k)*(  l.>dz)  +  dz*cl(  jz-f2,k) 
xr-xr-t-ar  ( k )  *p(  k) 
xi-xi-i-ai(k)*p(k) 

1  continua 

d2 5-cmplx ( xr , xi ) * f ac t 

return 

end 
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